#include "uhm.hxx"

#define N 1
#define DOF 3

using namespace uhm;

typedef ID_Node_H_<N>   Id_Nd;
typedef ID_Element_<N>  Id_El;
typedef LU_nopiv_Linal_ Uhm;

typedef Solver_<Id_Nd,Id_El,Uhm> Solver;

int main (int argc, char **argv) {

  initialize(argc, argv);

  // ** Environment set up
  Int_ 
    type = Mat_Base_::FLOAT,
    uplo = Mat_Base_::FULL,
    sym  = Mat_Base_::SYM,
    rhs  = 1,
    blk  = 256;

  Int_ 
    kind = Mesh_Base_::NODE_DEFAULT;

  // ** Time measuring begin
  time_begin();

  // ** Solver initialization
  Solver solver(type, uplo, sym, rhs, blk);

  typename Uhm::_A_ A[2];
  typename Uhm::_B_ B[2];

  for (Int_ i=0;i<2;++i) {
    A[i].set(type, uplo, sym, DOF*3, DOF*3);
    A[i].create_flat_without_buffer(); 
    A[i].create_flat_buffer();

    B[i].set(type, Mat_Base_::FULL, Mat_Base_::UNSYM, DOF*3, rhs);
    B[i].create_flat_without_buffer(); 
    B[i].create_flat_buffer();
  }

  random_matrix(A[0]);
  random_matrix(B[0]);

  if (get_verbose()) {
    std::cout << "A[0] = \n" << A[0];
    std::cout << "B[0] = \n" << B[0];
  }

  // ** ID objects
  Id_El el; el.set_id(0, 0);
  Id_Nd nd; nd.set_id(0, 0);

  // ** Supplying elemental matrices with mesh information
  std::vector< std::pair<Id_Nd,Int_> > nodes;
  nodes.reserve(3);

  time_in("solver:: Copying Matrices");
  solver.copy_in_begin();
  nd.set_id(0, 2); nd.set_dof(DOF); nodes.push_back(std::make_pair(nd, kind));
  nd.set_id(0, 1); nd.set_dof(DOF); nodes.push_back(std::make_pair(nd, kind));
  nd.set_id(0, 0); nd.set_dof(DOF); nodes.push_back(std::make_pair(nd, kind));

  solver.copy_in(el, nodes, 
                 A[0].flat().get_cs(), 
                 A[0].flat().get_buffer(), B[0].flat().get_buffer());
  nodes.clear();++el;

  nd.set_id(0, 4); nd.set_dof(DOF); nodes.push_back(std::make_pair(nd, kind));
  nd.set_id(0, 3); nd.set_dof(DOF); nodes.push_back(std::make_pair(nd, kind));
  nd.set_id(0, 2); nd.set_dof(DOF); nodes.push_back(std::make_pair(nd, kind));


  solver.copy_in(el, nodes, 
                 A[0].flat().get_cs(), 
                 A[0].flat().get_buffer(), B[0].flat().get_buffer());

  nodes.clear();++el;

  nd.set_id(0, 6); nd.set_dof(DOF); nodes.push_back(std::make_pair(nd, kind));
  nd.set_id(0, 5); nd.set_dof(DOF); nodes.push_back(std::make_pair(nd, kind));
  nd.set_id(0, 4); nd.set_dof(DOF); nodes.push_back(std::make_pair(nd, kind));

  solver.copy_in(el, nodes, 
                 A[0].flat().get_cs(), 
                 A[0].flat().get_buffer(), B[0].flat().get_buffer());
  nodes.clear();++el;

  solver.copy_in_end();
  time_out();


  if (get_verbose()) {
    std::cout << solver.work();
    std::cout << solver.mesh();
  }

  // ** Testing ops
  time_in("solver:: Mesh2Tree");
  UHM_ERROR((mesh2tree<Solver::_Mesh_,Solver::_Tree_>(solver.mesh(),solver.tree())),
            ">> mesh2tre Error");
  time_out();

  time_in("solver:: Separating Nodes");
  UHM_ERROR((op_leaf_to_me_par<Solver::_Tree_>(solver.tree().get_root(),
                                               &op_elt_separate_nodes<Solver::_Tree_>)),
            ">> Separating Nodes Error");
  time_out();

  time_in("solver:: Creating Matrices");
  solver.mesh().set_matrices(type, rhs, blk);
  solver.mesh().create_matrices();
  time_out();
  
  
  UHM_ERROR((op_leaf_to_me_seq<Solver>(&solver, solver.tree().get_root(), 
                                       &op_uhm_merge_A<Solver>)),
            ">> Simulating Merging Error");

  if (get_verbose()) {
    std::cout << "Mesh merged = \n" << solver.mesh();
  }

  // ** Randomize solution of 3rd element
  nd.set_id(0, 4); nd.set_dof(DOF); nodes.push_back(std::make_pair(nd, kind));
  nd.set_id(0, 5); nd.set_dof(DOF); nodes.push_back(std::make_pair(nd, kind));
  nd.set_id(0, 6); nd.set_dof(DOF); nodes.push_back(std::make_pair(nd, kind));

  solver.copy_out(nodes, B[1].flat().get_cs(), B[1].flat().get_buffer());
  nodes.clear();

  if (get_verbose()) {
    std::cout << "X copied  = \n" << B[1];
  }

  solver.clear();
  solver.flush();

  for (Int_ i=0;i<2;++i) {
    A[i].free_all();
    B[i].free_all();
  }

  time_end();

  finalize();

  return true;
}


